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Abstract: 

This paper describes the development of a nonlinear dynamic model for 
large oscillations of a robotic manipulator arm about a single joint. 
Optimization routines are formulated and implemented for the identification of 
electrical and physical parameters from dynamic data taken from an industrial 
robot arm. Special attention is given to the role of sensitivity in the 
formulation of robust models of this motion. The importance of actuator 
effects in the reduction of sensitivity is established and used to develop an 
electro-mechanical model of the manipulator system. 
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1 . Introduction 


The purpose of this research is to develop and investigate methods for 
identifying parameters in a dynamic model of a robotic manipulator. 
Identification routines of this type are important in the construction of 
control algorithms for manipulator systems [4], Because the parameter 
identification must be based on input and output data from an assembled 
manipulator, which acts under gravity and has possibly complicated joint 
friction, the dynamic model is a nonlinear differential equation, which must 
be solved numerically. 

The approach used to date is to employ a nonlinear search routine to 
minimize a quadratic fit-to-data criterion formed using the experimental data 
and the solution to the model equation. This method has been applied to a 
Unimation 600 Puma arm, with data obtained by F. W. Harrison in the 
Intelligent Systems Robotics Laboratory at the NASA Langley Research Center. 

Section 2 describes the mathematical model of the manipulator arm and the 
parameters to be identified. Section 3 describes the parameter identification 
scheme and the computer algorithms used. In Section 4, the experiment is 
discussed in more detail, along with some preliminary data reduction and 
analysis of the relationship between angular velocity and torque. Section 5 
presents an analysis of the sensitivity of the manipulator arm model to 
perturbations of uncertain parameters and initial conditions. We also discuss 
a method for reducing this sensitivity which in this application corresponds 
to the inclusion of a back electromotive force in the arm model. In Section 6 
we discuss the results of the parameter estimation routines for several models 


of robot arm friction 
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2. Manipulator Model 

To minimize the number of unknown parameters in each set of data, each 
experiment was performed with all but one manipulator joint locked. For each 
experiment, then, the model of the manipulator is a rigid arm that pivots 
about the one moving joint at a point 0. Thus, the arm m the model for a 
given experiment consists of several manipulator links, including the end 

effector, constrained to move as a rigid body. The equation of motion is 

(2.1) J0 - mgrsin 0 + f ( 0 ) = u(t) 

where 0 is the angle between the arm and the upward vertical and u is the 
control torque applied to the arm by the electric motor (actuator) at the 

joint in question. The damping term f(0) represents friction in both the 
joint and the motor; J is the moment of inertia about the appropriate joint, 

m is the mass of the arm, g is the acceleration of gravity, and r is the 

distance from 0 to the arm's center of mass. 

The angle between the arm and the vertical was measured at sampling times 
t^ and the sampling rate was 30 Hertz, so that 

(2.2) t = t... - t. = 1/30 sec. 

s i+i i 

We will denote this measured angle (i.e., the data) by y(t x ) to distinguish 
it from the solution to the model equation (2.1). 

The basic idea of the parameter identification scheme is to find 
parameters for (2.1) so that the solution to this differential equation 
matches the measured angle as closely as possible at the sampling times. 
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Because we cannot identify all of the parameters in (2.1) from the experiment 
described, we must define a minimal set of parameters for this model. 
Therefore, we rewrite (2.1) as 

(2.3) 9 - a sin 0 + f(c^, c^, 9) = 0u(t) 

where a = mgr/J, 0 = 1/J, and we have parameterized the the damping term 
f ( 9 ) in (2.1) as f(c^,C2,0). The damping term f ( 0 ) may include various 
forms of dissipation: linear damping, nonlinear damping, and Coulomb 

friction. Our best results have come with (sometimes piecewise) linear damping 
and quadratic damping in combination with a linear damping term resulting from 
back electromotive force. 

We will refer to the set of parameters in (2.3) by the parameter vector 

(2.4) q = [ct 0 c^ c 2 • • •]. 

3. Parameter Identification 

An experiment performed on a time interval [tg,t^] yields data u(t^) 
and y(t i ), t^ = tg , t Q + t g ,»**, t f . With the known command torque u(t) 
and a set of trial parameters, we solve (2.3) on the interval [tg,t^] and 
form the fit-to-data criterion 

(3*U J(q) = l [6(t i ) - y(t i )] 2 . 

The parameter identification then consists of finding the parameter vector q 
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to minimize J(q). Usually, we take the initial time t^ > 2 sec. because we 
suspect some error in the data near the beginning of the experiment due to 
transients in electronics. Therefore, in some cases we know that the initial 
angular velocity is zero, but in most cases we must estimate it using finite 
differences obtained from the position measurement. 

To solve (2.3), we use a fourth-order Runge-Kutta algorithm with variable 
step size [3]. We tried using the numerical integrators DGEAR and DVERK in 
the IMSL library, but both of these routines often hung up — i.e., the step 
size was reduced to zero — where the manipulator arm turned. This was 
especially troublesome for models with piecewise continuous damping and 
Coulomb friction. The step-size control in our final Runge-Kutta routine does 
not allow the step size to fall below a specified minimum. 

For minimizing J(q) we used the subroutine ZXSSQ from the IMSL library, 
which is a Levenberg-Marquardt algorithm [2] that approximates gradients by 
finite differences. It also estimates the Hessian. Hence we assume certain 
smoothness and local convexity of J(q) and the performance of the algorithm 
indicates that these assumptions are valid. 


4. Data Collection and Analysis 

Experimental data was collected by F. W. Harrison in the Intelligent 
Systems Robotics Laboratory at NASA Langley Research Center. The subject of 
the experiments was a UNIMATE PUMA industrial robot with six degrees of 
freedom. A schematic [1] of the robot arm with rotational joints is shown in 
Figure 1. The experiments described below were performed by rotating only the 
shoulder (joint 2) with all other joints locked in a collinear position. 
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a. Static Experiment 

The purpose of this experiment was to establish a relationship between the 
current delivered to joint 2 and the torque exerted on the robot arm. The arm 
was placed in a horizontal position and the force exerted by the arm at a 
fixed distance from the joint was measured at a sequence of motor current 
levels. The results of this experiment are shown in Figure 2. Because of the 
linear relationship of these quantities, it was decided to use the measured 
motor current data as the torque input in our dynamic models of the arm 
motion. 

b. Dynamic Experiment 

The purpose of this experiment was to gather input and output data for the 
dynamic model described in Section 2. The arm was initialized in a vertical, 
upright position and then commanded to rotate about joint 2 through an angle 
of approximately 90 degrees in both directions. During this oscillation, 512 
measurements of the joint angle in radians and the motor current as measured 
by the voltage drop across a known resistance were taken at a frequency of 30 
Hertz. These input and output data are illustrated m Figures 3 and 4, 
respectively. The angular velocity of the robot arm, calculated by backward 
differences, is shown in Figure 5. 

c. Command Torque Synthesis 

This data analysis is designed to recover the square-wave commanded 
voltage across the motor terminals. If only back electromotive force is 
included in the motor model, then this commanded voltage is the sum of the 
voltage drop across the motor resistance (Figure 3) and the angular velocity 
(Figure 5) multiplied by the back emf constant. Figure 6 shows the results of 
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this computation for a back emf constant of 1.5, estimated by trial-and-error . 
In Section 5 a square-wave approximation of Figure 6 is used as input to an 
alternate model of the robotic system which includes back emf effects. 


5, Sensitivity Analysis 

In this section we discuss the sensitivity of the solution of the 
nonlinear model 


{ 0 - a sin 0 + f(c^, C 2 » 9) = Bu(t) 

9(0) = 0 Q 

0(0) = o)q 

with respect to small perturbations of the initial velocity Wq and the 
unknown friction parameter c^ . 

In some applications of parameter estimation, moderately high parameter 
sensitivity is advantageous in that it allows the unknown parameters to be 
estimated with a greater degree of certainty for a given level of noise in the 
output data. However, in parameter estimation for simulation, a sensitive 
model will yield poor simulations when the unknown parameters are subjected to 
slight variations due to modeling errors or external factors. For this reason 
it is preferable to have a mathematical model which is relatively stable with 
respect to perturbations of the parameters. 

The same reasoning applies to perturbations of initial conditions. In 
this application only the initial position can be measured directly for a 
selected subinterval of the data. The angular velocity, computed by a finite 
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difference, is subject to considerable error. Therefore, it is desirable to 
have a simulation model which is not subject to high sensitivity with respect 
to this initial condition. 

Unfortunately, numerical testing indicates a substantial degree of 
sensitivity on both parameters and initial conditions for the model (5.1) with 
measured motor current as input. Figure 7 indicates the effect on the output 
of a small perturbation of the friction coefficient cp and Figure 8 shows a 
similar comparison for a small perturbation of the initial velocity Wq. In 
fact this figure indicates an almost chaotic dependence of this model on the 
initial angular velocity. 

We therefore undertake a mathematical method for reducing this sensitivity 
which in this application has physical implications as well. The basic idea 
is to increase the damping in the system by adding a term of the form 0k0 to 
both sides of the differential equation. Setting v(t) = u(t) + k0(t) yields 

(5.2) 0 - a sm 0 + [f(cj,C2,0) + 8k0 ] = f3v(t) 

which is of the same form as (5.1) except that the damping term has been 
increased at the expense of changing the input function. If, in the parameter 
estimation procedure, the new input v(t) can be either measured directly or 
synthesized from data, then equation (5.2) is a model which may possess 
greater stability with respect to initial conditions and parameter values. 
This factor tends to yield numerical solutions for the state which are more 
reliable over long time intervals and therefore lead to more robust parameter 


estimates . 
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A numerical test of this procedure is shown in Figures 9 and 10. These 
figures show the output of a model based on (5.2). The input v is the 
synthesized commanded motor terminal voltage (Figure 6). The damping 
parameters are greater because they include back emf effects. The numerical 
results show significantly lower sensitivity on friction parameters (Figure 9) 
and initial angular velocity (Figure 10). 

The physical implications of the mathematical procedure for this system 
have already been indicated. By using the commanded motor terminal voltage as 
input rather than the torque delivered to the joint, one arrives at a more 
stable model of the robot dynamics. In effect, the model which includes back 
emf takes advantage of a natural damping in the electro-mechanical system. 
One can obtain good fit-to-data over short time intervals for the mechanical 
system alone, but the stability provided by this effect is lost. 


6. Numerical Results 

The iterative parameter estimation routine described in Section 3 was 
applied to the model (5.1) on the time interval [1.67, 15.0]. The results are 
given in Table 1. Three alternative friction models were employed: linear 

friction given by 

(6.1) f(x) = cx, 

quadratic friction of the form 


( 6 . 2 ) 


f (x) = c^x + C 2 x| x| , 
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and piecewise linear, direction-dependent friction model of the form 


(6.3) 


f(x) 


I ^x, x > 0, 

c 2 x, x < 0. 


As shown in Table 1, the cost function (3.1) is quite large for all three 
friction models and the iterative method does not converge. 

The results of the same procedure for the desensitized model (3.2) with 
input given in Figure 6 are shown in Table 2. The friction values are larger 
because they include the effect of back eraf. The procedure shows convergence 
to relatively low cost values for each of the three friction models. The 
direction-dependent friction model (6.3) shows the most rapid convergence to 
the lowest cost value. While the quadratic friction model (6.2) eventually 
obtains a low cost, it converges much more slowly and alters the physical 
parameters a and ft to values which are quite different from those obtained 
by models (6.1) and (6.3). The solid graph in Figure 10 show the output of 
the most accurate model in Table 2. This graph is almost indistinguishable 
over this time interval from the graph of the measured robot arm motion 
(Figure 4). 


7. Conclusion 

Our experience indicates the importance of actuator effects in the 
development of robust dynamic models for the motion of a robotic manipulator 
arm. The inclusion of natural damping due to back eraf effects improved the 
performance of both the numerical integrator for solving the nonlinear 
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equation of motion and the numerical optimizer for estimating parameters. 
Among the friction models we studied, the model allowing direction-dependent 
damping coefficients was the most successful. In continuing research, we plan 
to combine this model with a first-order dynamic model of the actuator to 
study more complex motions of the robotic manipulator arm. 
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Table 1. Numerical Results for Model (5.1). 


Linear friction model (6.1) 


iterate 

a 

8 

c 

cost 

0 

13.39 

21.00 

1.17 

82300. 

1 

16.33 

24.07 

1.08 

65600. 

2 

15.51 

22.64 

1.12 

137000. 

3 

11.72 

13.64 

1.18 

6110. 

4 

13.18 

11.58 

3.23 

5410. 

5 

14.14 

10.31 

4.03 

5240. 


Quadratic friction model 

(6.2) 




iterate 

a 

6 

C 1 

c 2 

cost 

0 

13.39 

21.00 

1.17 

0.00 

82300. 

1 

31.55 

38.38 

3.13 

0.55 

152000. 

2 

26.72 

23.31 

5.47 

0.58 

103000. 

3 

180.4 

130.5 

27.2 

2.96 

56800. 


Piecewise linear 

friction 

model (6.3) 




iterate 

a 

8 

C 1 

c 2 

cost 

0 

12.05 

21.42 

2.59 

4.87 

450. 

1 

11.23 

20.73 

3.49 

4.67 

848. 

2 

10.65 

21.26 

4.22 

5.50 

720. 

3 

10.94 

21.40 

4.08 

8.79 

734. 

4 

10.40 

20.98 

4.09 

1.86 

25000. 

5 

11.28 

20.66 

3.64 

4.07 

15300. 

6 

12.26 

19.70 

1.65 

6.52 

14400. 
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Table 2. Numerical Results for Model (5.2). 


Linear friction 

model (6.1) 




iterate 

a 

6 

c 

cost 

0 

12.00 

20.00 

35.00 

301. 

1 

10.91 

20.99 

34.04 

25.2 

2 

11.09 

21.44 

33.55 

5.14 

3 

13.25 

21.41 

33.42 

2.85 

4 

13.57 

21.44 

33.34 

2.60 

5 

13.39 

21.01 

32.67 

2.59 


Quadratic friction model 

(6.2) 




iterate 

a 

8 

C 1 

c 2 

cost 

0 

12.00 

20.00 

35.00 

0.00 

301. 

1 

11.01 

20.85 

34.19 

-0.72 

14.8 

4 

14.19 

21.93 

32.87 

1.08 

2.26 

7 

18.74 

25.20 

27.19 

10.59 

0.89 

14 

20.23 

26.17 

24.07 

14.81 

0.66 


Piecewise linear 

friction 

model (6.3) 




iterate 

a 

6 

C 1 

c 2 

cost 

0 

12.00 

20.00 

35.00 

35.00 

301. 

1 

11.27 

20.70 

34.30 

35.70 

5.00 

2 

11.31 

20.91 

34.20 

35.68 

0.21 

3 

11.81 

21.09 

34.21 

35.49 

0.116 

4 

12.05 

21.42 

34.72 

36.00 

0.112 
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Robot arm with rotational joints [1] 
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Motor Current 


Figure 2. Static torque regression. 


Correlation coefficient is 0.999. 



Motor Current (volts) 
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Time (sec) 


Figure 3. Input. The motor current was measured by the voltage drop across a 


known resistance. 



Angular Position (radians) 
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Time (sec) 


Figure 4. Output. 



Angular Velocity (rad/sec) 
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Time (sec) 


Figure 5. Differentiated output. 



Synthesized Motor Terminal Voltage 
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Figure 6. Synthesized input. This graph is the sum of Figure 3 and 1.5 times 


Figure 5. 
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Time (sec) 

Figure 7. Parameter sensitivity for model (5.1) with friction model (6.3). 

Common parameters are a = 12.74, 3 = 21.209, and c 2 = 4.301. The 
perturbed parameter is c L with Cl = 2.296 in the solid graph and 
c^ = 2.306 in the dotted graph. 




Angular Position (radians) 
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Time (sec) 

Figure 8. Initial angular velocity sensitivity for model (5.1) with friction 
model (6.3). Parameters are a = 12.74, g = 21.209, c^ = 2.296, 
and C£ = 4.301. The angular velocity is = 0.63 in the solid 
graph and o)q = 0.64 in the dotted graph. 



Angular Position (radians) 
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Time (sec) 

Figure 9. Parameter sensitivity for model (5.2) with friction model (6.3). 

Common parameters are a = 12.05, 0 = 21.42, and c 2 = 36.00. The 

perturbed parameter is c x with c L = 34.72 in the solid graph and 

c 2 = 33.72 in the dotted graph. 




Angular Position (radians) 
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Time (sec) 

Figure 10. Initial angular velocity sensitivity for model (5.2) with friction 
model (6.3). Parameters are a = 12.05, g = 21.42, and 

= 34.72, and C 2 = 36.00. The initial angular velocity is 
<i)q = 0.76 m the solid graph and oJq = 1.76 in the dotted graph. 
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